!************************************************************
! this subrouitne calculate invese Symmetry point of
! isentropic supersonic flow
!************************************************************
! Created by : B. G.
! Date       : 2015-08-13
! Revised    :
!************************************************************
! Nomenclature: 
!   Vm~Mm   properties along r-c-line 14
!   p00~a00 properties along streamline 34
!   lm~T02  coefficients of functions, *m--r-c-line, *0(*)-streamline
!   *4f     properties of previous loop of point 4
!   i       loop indicator
!   
!************************************************************
! Warning :
!   1. angle in rad, not degree
!************************************************************

subroutine SymmetryP()
   !      1
   !
   !  3      4()
   use VariableDef
   implicit none
   ! given variables
   ! real*8::x1,y1,V1,theta1,p1,rho1,x3,y3,V3,theta3,p3,rho3
   ! temporary variables
   real*8::Vm,pm,rhom,thetam,ym,alpham,Mm
   real*8::p00,rho0,V0,a00
   real*8::lm,Qm,Sm,Tm,R0,A0,T01,T02
   real*8::x4f,p4f,V4f,rho4f
   integer::i
   ! calculate variables
   ! real*8::x4,C4,p4,rho4,T4,M4,Vv4,Vu4
   
   p4=p3
   V4=V3
   rho4=rho3
   y4=0.0
   theta4=0.0
   
   Vm=V1
   pm=p1
   rhom=rho1
   thetam=theta1
   ym=y1
   
   x4f=0.0
   p4f=0.0
   rho4f=0.0
   V4f=0.0
   
   do i=1,icor
      ! determine coefficients along r-c-line
      call Thermo(pm,rhom,Vm,a00,T4,Mm)
      alpham=asin(1/Mm)
      lm=tan(thetam-alpham)
      Qm=sqrt(Mm**2-1)/(rhom*Vm**2)
      Sm=delta*sin(thetam)/(ym*Mm*cos(thetam-alpham))
      ! determine x4 and Tm
      x4=x1-y1/lm
      Tm=-Sm*(x4-x1)+Qm*p1-theta1
      ! determine coefficients along streamline
      p00=(p3+p4)/2
      rho0=(rho3+rho4)/2
      V0=(V3+V4)/2
      call Thermo(p00,rho0,V0,a00,T4,M4)
      R0=rho0*V0
      A0=a00**2
      T01=R0*V3+p3
      T02=p3-A0*rho3
      ! determine p4, rho4, V4
      p4=(Tm+theta4)/Qm
      V4=(-p4+T01)/R0
      rho4=(p4-T02)/A0
      ! justify if point 4 arrived
      if ((abs(x4-x4f) <= el) .and. (abs(p4-p4f) <= ep*p4f) .and. &
         &(abs(rho4-rho4f) <= erho*rho4f) .and. (abs(V4-V4f) <= eV*V4f)) exit
      
      pm=(p1+p4)/2
      rhom=(rho1+rho4)/2
      Vm=(V1+V4)/2      
      thetam=(theta1+theta4)/2
      ym=(y1+y4)/2
   
      x4f=x4
      p4f=p4
      rho4f=rho4
      V4f=V4
   end do
   call Thermo(p4,rho4,V4,a00,T4,M4)
   Vv4=V4*sin(theta4)
   Vu4=V4*cos(theta4)
   
end subroutine SymmetryP
